High-temperature surface superconductivity in rhombohedral graphite 
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Surface superconductivity in rhombohedral graphite is a robust phenomenon which can exist even 
when higher order hoppings between the layers lift the topological protection of the surface flat band 
and introduce a quadratic dispersion of electrons with a heavy effective mass. We show that for 
weak pairing interaction, the flat band character of the surface superconductivity transforms into 
a BCS-like relation with high critical temperature characterized by a higher coupling constant due 
to a much larger density of states than in the bulk. Our results offer an explanation for the recent 
findings of graphite superconductivity with an unusually high transition temperature. 

PACS numbers: 74.20.-z, 74.20.Pq, 74.70.Wz 



A low critical temperature of conventional supercon- 
ductors results from a constant density of states (DOS) 
due to a linear-in-momentum electronic spectrum near 
the Fermi energy. With a higher-order dispersion, the 
relation between the critical temperature and the cou- 
pling constant becomes stronger, boosting the supercon- 
ductivity. The extreme case would be a completely dis- 
persionless energy spectrum, a fiat band, which has been 
predicted in many condensed matter systems, see e.g. 
Refs. [ll-H. In some cases the flat bands are protected 
by topology in momentum space; they emerge in gap- 
less topological matter [5 -13 1 . A singular DOS associated 
with the dispersionless spectrum was recently shown 
to essentially enhance the transition temperature open- 
ing a new route to room-temperature superconductivity. 

The problem is to find the metal with such a higher- 
order dispersion around the Fermi sea. Refs. and 14 
have shown that within the nearest-neighbour approxi- 
mation, rhombohedral graphite (RHG) has topologically 
protected surface states with a flat band at the Fermi en- 
ergy, and these surface states support high-temperature 
superconductivity where the superconducting order pa- 
rameter is concentrated around the surfaces. A flat band 
forms out of a low dispersive band that appears on the 
surface of a multilayered rhombohedral graphene struc- 
ture with a large number of layers. The corresponding 
critical temperature depends linearly on the pairing in- 
teraction strength and can be thus considerably higher 
than the usual exponentially small critical temperature 
in the bulk. Flat-band superconductors can carry quite 
high surface supercurrent with the critical value propor- 
tional to the critical temperature |15|. 

Experimental evidence of a high-temperature super- 
conductivity in graphite in the form of a small Meiss- 
ner effect and of a sharp drop in resistance appeared in 
the literature during past years [l^ 12 1- Recently, these 



cations of even room-temperature superconductivity in 
specially prepared graphite samples |19j . The enhanced 
superconducting density has been also reported on twin 
boundaries in Ba(Fei_^Co2;)2As2 [l^]. In this Letter, 
we argue that the high-temperature superconductivity in 
graphite can be related to surface superconductivity that 
may form either on the outer surfaces of the sample or 
on twin boundaries of or on grain boundaries between in- 
clusions of RHG. We demonstrate that the surface super- 
conductivity is a robust phenomenon which survives even 
when the topological protection of the flat band itself is 
lifted. In particular, the next-nearest neighbour hoppings 
in RHG can break the exact topological protection and, 
therefore, the flat-band mechanism of superconductivity 
could be destroyed. We study these higher-order interac- 
tions 2l| and show that, though breaking the fiat-band 



scenario at sufficiently low values of the coupling energy, 
they provide another mechanism of surface superconduc- 
tivity which is of the BCS type but still has a much larger 
coupling constant than the usual superconductivity in 
bulk, thus favoring high-temperature superconductivity. 
The enhanced coupling constant comes from a high DOS 
associated with a heavy effective mass of surface quasi- 
particles emerging on the background of the pre-existing 
fiat band. Our results help to identify the regime of pa- 
rameters where extremely high-temperature surface su- 
perconductivity may be found. 

Electron dispersion in RHG. The RHG lattice and 
the tight-binding couplings are depicted in Fig. [1] We 
label the layers (starting from the bottom) by n, the 
atoms A in layer n are on top of atoms B in layer n — 1 , 
and the vector between the layers is d. In our estimates, 
we use the tight-binding parameters denoted in Fig. [1] 
They satisfy 70 3> 71 



73 > 74 [22| 



findings have been ratified by observations of zero resis- 
tance in graphitic samples up to 175 K [isj and indi- 



In the tight- 
binding numerics below, we use the values 70 = 2.58 eV, 
71 = 0.34 eV, 73 = 0.17 eV, and 74 = 0.04 eV, which 
give the best fit to the density functional theory (DFT) 
calculation of the surface state dispersion (see Appendix). 
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FIG. 1: Rhombohedral graphite. The black and gray atoms 
correspond to A and B sites, respectively. 



The RHG is a multilayered graphene structure. The 
conical spectrum near the Dirac points K and K' of the 
Brillouin zone of a single-layer graphene (for details, see 
review 



22l | and references therein) is transformed into 



low-dispersion, low-energy bands (sec Fig. ^ which de- 
termine the unique features of this system. Since we are 
interested in low energies we concentrate on the in-plane 
momenta p = [px, Py) close to one of these Dirac cor- 
ners. A standard Fourier series expansion near K yields 

H 



N 



Hk^Y^Y. il(P)ffm„(K,p)Vi„(p), (1) 



p 7n,n—l 



where i/™„(K, p) = J2t=a -ffmn(K, p) and 
H^l{K,p) = -71 [e"'^a-+(5„,„+i -He*-&(T_(S,„,„_i] 

H^l{K,p) =74fF [e'^ P-Sm.n^l + e~''-ip+Sjn,n+l] ■ 

Here 73 = 73/70, 74 = 74/70, P± ^ Px ± ipy = pe='=*'^, 
and I'F = 3ao7o/2?i. The Pauli matrices & and 2ct± = 
ax ± itTy act on pseudo-spinors tpn = {ipn , V'n)"^, V'i = 
(Vi* , Vr), where ^^^ = ,p^,^l= e^/^^- 

To construct the associated Bogoliubov-de Gennes 
(BdG) Hamiltonian for the superconducting state we 
need also the time-reversed "hole" Hamiltonian for the 
Dirac point K. It follows from the particle Hamiltonian 
in a vicinity of the opposite Dirac point — K which is 
equivalent to K'. The wave function V-'k of ^ hole excita- 
tion near K is ip^ = i/'Ik- ^Jne can check that the hole 
Hamiltonian is i7^„(K,p) = i/*„„(— K, — p). Therefore, 

N 

^K=E E ^^')t(p)^mn(K,p)Vii"nP)- (2) 
p m.n — 1 

In what follows, we denote the electron wave function by 
Un = ipn and the hole wave function by u„ = ipn- ^^c 



FIG. 2: Cuts of the 3d spectrum for p < pfb = 'yi/vp along 
the K-F (K-M) on negative (positive) x-axis; A'' = 5 (red), 
iV = 10 (blue), and A'' = 20 (black) graphene layers. The 
DFT calculations are shown in full, tight-binding in dash- 
dotted and Eq. ((Tj) in dashed lines (the latter plotted up to 
their regimes of applicability). The inset is a zoom- up of the 
low-energy region with tight-binding (dash-dotted) and ana- 
lytical (dashed) curves. The deviations between the dashed 
and other lines show up when becomes dominant in Eq. ((7|, 
and are partially due 73 neglected there. 



next section we study the normal state using Hamilto- 
nian Eq. ([l} and the relative magnitudes of the coupling 
constants listed above. The results of the numerical solu- 
tion are displayed in Fig. [21 along with the corresponding 
analytical approximations and DFT calculations. 

Low-energy spectrum in the normal state. The 
Schrodinger equation takes the form 



E ^nm(K, p)Um(p) = (e + Ai)""(p) 



(3) 



The energy is measured from the chemical potential /i. 
The energy spectrum in bulk is obtained by ignoring the 
outermost layers n = 1 and A'^ and using the ansatz 
Un e'^^"*", where p^ is out-of-plane momentum. For 
zero doping ^ = 0, the Fermi surface is determined by 
e{p,pz, <p) = 0. If 73 = 74 = 0, the Fermi surface shrinks 
to a spiral line vfP = li t 4> = Pzd + 7r/6. Projection of 
this spiral onto the momentum plane g = determines 
the area of a flat band for surface states Q in the limit 
N ^ 00. If only 74 = while 73 7^ 0, equation e = 
for /i = can still be shown to give a Fermi surface in 
the form of a (corrugated) spiral whose projection de- 
termines the surface flat band 2l|. This is because the 



interaction comes with matrices Ux and ay so that the 
full Hamiltonian obeys the same anti-commutation rule 
[ct^, {H'^^^ + H^^^ + H^^^)]+ = as the initial Hamiltonian 
Hio) +hW, This preserves the same topological invari- 
ant and the same topology of the Fermi surface 0] . Since 
73 docs not affect the presence of a flat band, we restrict 
our analytical consideration to the case when only 74 is 
nonzero while 73 = for simplicity. However, our numer- 
ical analysis is carried out using the full Hamiltonian. 
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Surface states have complex Pz = p'z ~^ w'z ^'^'^ decay 
into the bulk. For 73 = Eq. ([3]) in the particle channel, 

+74 [e'^WFP-Un+l + e'^'^VFP+Un-l] = (e + /i)'Un , (4) 

for low energies has a solution in the form 



1 



c 



(5) 



where p = p/ppB, Pfb ~ ^i/vp and 

C = p[{e + m)/7i - 74(p' + l)]/(p' - 1) . (6) 



Here the out-of-planc momentum p'^d 



tt/6 while 



g±Pj <^ = The overall normalization C is found from 
dX]^=i[Tr^JjU„] = 1. For large N this gives |Cp = 
provided \A+\'^ + [^-p = 1. 
At the outermost layers, the terms with ipQ and tpN+i 
in Eq. (|4]) disappear. The components which do not 
have 7i couple the constants and A~ in Eq. ([5]) and 
determine the energy of the surface states 



ep= fip± (1 



f) 



(7) 



for ^p, e <C 71. Here 

fip =■ p^ /2m* - ^ , m* = 71/ (4741;^) . (8) 

The interaction 74 breaks the symmetry between the con- 
duction and valence bands in a way similar to a shift in 
/i due to doping. The spectrum Cp has a quadratic dis- 
persion with the effective mass m* on a background of a 
much weaker high-order dispersion ^p. The latter trans- 
forms into a flat band ^p = with a radius p < pfb for 
an infinite number of layers, N 00. The effective mass 
is much larger than the characteristic band mass in 
3D graphite. Indeed, we have m* /m^ ~ 71/74 where we 
estimate / (m^ag) 70 as the conduction band width 
in graphite. We see that m* /m^ ^ 1. This dispersion is 
compared with the results of numerical diagonalization 
of H{K, p) in Fig. fusing 73 = 0.066, and 74 = 0.016. 

BdG equations are constructed using the particle and 
hole Hamiltonians ^ and ([2]) coupled through the super- 
conducting order-parameter field A. As distinct from the 
quasiparticle energy measured from the chemical poten- 
tial upwards, E = p + e, the energy of holes is measured 
from fi downwards, E = p — e. We have 



E 



T3 



n 



. (9) 



Here we introduce objects in the Nambu space 



" " I A* ' ' " " 



T3 = 



1 

-1 



Each component of the Nambu vector ^„ is a pseudo- 
spinor. For analytical consideration we assume that 
A„ = for n ^ 1,A^. This is justified by the numeri- 
cal solution of the self-consistency equation using the full 
BdG equations [lj,l2l|. In this case, Eq. © for n 1 , 
does not contain A, so that one can use the normal- 
state solution, Eq. ([5]), where we have Nambu vectors 
A^ = {A^, B^)^ instead of the corresponding scalars, 
and the Nambu matrix C, Eq. ([6]), with f^e instead of e. 

At the outermost layers, the terms with ug, and 
iiM+i, vpf+i in Eq. ([9|) disappear. The components which 
do not contain 71 yield 



= (? - f3Pp)A+ - AiA+ , 
f3^pA+ = (? - f3flp)A- - An A- , 



(10) 
(11) 



where e = e (l — p^) ^, jlp = fip (l — p^) ^ ■ Equations 



([TUl) . ([TT1) provide the surface-state spectrum and deter- 
mine four independent surface states. 

If Ai = Ajv, the spectrum is = (/ip ± ^pf + |Ap. If 
the number of layers A^ is large, ^p — for p < pfb, the 
two surface states decouple 



\A\l 



\A\% 



In this case, Eq. ^TU\i at layer n ~ 1 yields 
B+ or A+ ^V, B+ where 



U 



[l+~Pp/e]^ , V 



[1 - 



(12) 

= u, 

(13) 



Surface superconductivity. The surface states dis- 
cussed above form a basis for the superconducting gap 
localized near outer surfaces. Within the mean-field ap- 
proximation, the gap at layer n is determined by the 
self-consistency equation. As was shown in Ref. [iJ, the 
surface states dominate due to a much larger DOS. For a 
large number of layers when ^p = 0, the self-consistency 
equation for the gap at the surface takes the form 



FB 



^ii^tanh^ 



(14) 



Here we used Eq. ((13)) to find the integrand, e is one of the 
spectral branches in Eq. p2p . the integration is carried 
out over momenta within the fiat band, p < PfBj and W 
is the 3D coupling potential. This is the central result 
of our Letter. The superconducting coupling is described 
by the energy g — (W / d)p']^ g / ti^ . It can also be expressed 
in terms of the usual BCS coupling constant A = v^W 
where V3 ~ m^p^p /2Tr^h^ is the 3D density of states and 
P3P is the Fermi momentum in 3D graphite. Assuming 
the conduction band width in 3D graphite of the order 
of 70 we have 3/71 A(7i/7o) if h/aoPsF ^ 1- 

The overall behavior of A vs. the coupling energy is 
plotted in Fig. (3] Let us consider first the resulting A 
for zero doping /.j = 0. The quadratic dispersion due 
to 74 comes with an energy scale a = 27471, which de- 
termines a crossover between exponentially suppressed 
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FIG. 3: Self-consistent surface gap vs. coupling g. The black 
line (in the middle) shows the results based on the exact diag- 
onalization of the BdG equations for A'' = 20 with jJ. — 0, the 
blue (bottom) line is Eq. H14[) at /i = and the red (top) line 
corresponds to fi = /lopt that maximizes A for given g. For 
large g, the gap tends towards the flat-band limit A oc g. For 
g < 4na (upper inset), the gap is exponentially suppressed, 
A (X exp{—4:TTa/g). The lower inset shows the (normalized) 
gap as a function of fi for a few values of g. 



and flat band superconductivity. For g ^ Ana and for 
zero doping, Eq. yields the flat- band result [l^ . 

A = g/Sir ioi T — 0, and the critical temperature satis- 
fying A = iksTc- Due to its linear dependence on the 
interaction strength, the critical temperature is propor- 
tional to the area of the flat band and can be essentially 
higher than that in the bulk. Doping in the flat band 
regime destroys the surface superconductivity . Both 
Ao and Tc vanish at the critical doping level |^| = 2kBTc. 

For g <C 47ra the weak dispersion Eq. ([7]) with a heavy 
mass TO* dominates. The integral in Eq. (jl4p is logarith- 
mic which results in a BCS-like expression (for T = 0) 



A = [a^/ia - M)]e"^^^' , A2 = g{l - ti/af/Ana 

where a = 27174. The estimate for g gives A2 ~ 
^(71/74)- This is a much larger coupling constant than X 
for bulk superconductivity. The gap disappears at ^ = a. 
The crossover from the BCS-like to the flat-band regime 
occurs at g Ana, and A ^ a. The coherence length 
^0 = hvg/A ~ ao (70/71)6^^'^^ is much longer than the 
interatomic distance oq. These analytical results arc 
compared in Fig. [3] to the numerical solution of the self- 
consistency equation using the full Hamiltonian, Eq. ^ . 
In contrast to the flat-band and the BCS regimes, the 
gap in the intermediate region a ~ g is enhanced by an 
optimum doping, i.e., the critical temperature is very sen- 
sitive to the presence of impurities (lower inset in Fig. [3]) . 
This complies with the reports of high-temperature su- 
perconductivity in doped graphite, Ref. 

Effect of fluctuations. The quality of the mean-field 
approximation used above is determined by the Ginzburg 
number Gi which is a measure of the relative magnitude 



of order-parameter fluctuations. For usual 3D supercon- 
ductors Gi ^ 1 due to a small ratio of the critical tem- 
perature to characteristic energy of electrons (i.e., the 
Fermi energy). Here we demonstrate that the mean field 
approach also works well when the quadratic dispersion 
dominates over the flat band. In this case the fluctu- 
ation free energy density for T not too close to Tc is 
Fi - U2AI/2, where 1/2 = m* /2nh^ is the 2D DOS and 
the effective mass to* is determined by Eq. jS]). The en- 
ergy of an area tt^q with a radius of ^0 = ^"'^'ff^cT^ is ~ 
TT^gFi = 747i(AJ/Aq), where Aq is the mean-field gap. 
Since J"i - T we find A?/A^ = Gi - Tc/7471. When the 
quadratic dispersion dominates, one has T^ <S 7i74 with 
the Ginzburg number Gi = e~^/^^ ^ 1, thus the average 
fluctuation of the order parameter is small compared to 
its mean-field value. However, at the crossover to the flat 
band regime, the fluctuation becomes of the same order 
as the mean-field Aq ^ 7174. Therefore, the mean-field 
approach is not exact for the flat-band regime. Never- 
theless, it is used here as an initial step towards a full 
theory of high temperature surface superconductivity. 

Summary. Rhombohcdral graphite is a promising 
candidate for high-temperature surface superconductiv- 
ity due to its (approximate) topologically protected flat 
band. Besides surfaces, similar type of superconductivity 
may arise around stacking faults and interfaces between 
differently stacked regions of graphite, as long as the sys - 
tem contains more than a few layers of RHG regions [21 1 . 
Recent observations of high-temperature superconductiv- 
ity in graphite [l6l - lil9i | are compatible with surface or in- 
terface superconductivity described by our theory if there 
are RHG regions embedded inside otherwise Bernally 
stacked graphite. Our predictions can be used for search 
or for an artiflcial fabrication of layered and/or twinned 
systems with high- and even room-tcmperature super- 
conductivity. With the hopping parameters used above, 
the crossover between the flat-band and BCS-like regimes 
takes place around gc ^ Ana « 0.397i « 0.15 eV (see 
Fig. [3] and Appendix) corresponding to the mean-fleld 
Tc{gc) ~ 20 K. For g > g^ we thus find T^ - (g/71) x 50 
K. This is much greater than the expected gap in the 
bulk for the same magnitude of coupling. 
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APPENDIX: COMPUTATIONAL DETAILS 

The density-functional theory calculations on rhom- 
bohedral graphene slabs were performed using the all- 
electron FHI-aims code We used the Perdew-Burke- 
Ernzerhof (PBE) exchange-correlation functional for all 
calculations, and took into account the van der Waals in- 
teraction using the approach by Tkatchenko and Scheer 
0. The Brillouin zone was sampled using a 48x48x1 
/c-point grid. "Tight" basis defaults as defined in the 
FHI-aims distribution were used for the all-electron de- 
scription of the carbon atoms. 



The in-plane lattice parameter was optimized for 
monolayer graphene and it was found to be 2.466 A. The 
same lattice parameter was used for the rhombohedral 
slabs but the distance between adjacent layers was al- 
lowed to relax freely until forces acting on atoms were 
less than 0.001 eV/A. This resulted in variations of the 
order 0.01 A in the interlayer distances, the outer lay- 
ers being slightly expanded from the optimized interlayer 
distance 3.332 A of rhombohedral graphite bulk. The use 
of the lattice parameter of graphene leads to negligible 
errors, as the difference to the optimized bulk parameter, 
2.462 A, is very small, and the multilayer structures are 
expected to interpolate between the bulk and monolayer 
limits. 



FITTING THE TIGHT-BINDING MODEL 

The band structures were calculated along the T~K — 
M direction in the Brillouin zone, sampling the lines be- 
tween (0.3,0.3) and (1/3,1/3), as well as between (1/3 
and 1/3) and (0.3,0.35) using 200-^200 fc-points. The 
band structure in the vicinity of the if-point was fitted to 
a tight-binding model, in which hoppings between near- 
est neighbors in-plane (70) and out-of-plane (71), as well 
as between next-nearest neighbors out-of-plane (73 and 
74) were included. The fit was concentrated to the region 
around K using a Gaussian function to weigh the squared 
error between the tight-binding and DFT bands, and the 
minimal energy of the parabolic region was shifted to 
match the corresponding DFT energy. The width of the 
Gaussian was chosen such that the weight outside the flat 
band region was practically zero. Fig. Hl^a) demonstrates 
a comparison between the DFT band structure and the 
corresponding tight-binding fit, as well as the shape of 
the weighing function on a 20-layer slab. 

Two parameter sets with comparable agreement 
around the iiT-point were found when the width of the 
weighing Gaussian was altered. As Fig. SJ^b) and (c) 
show, one of them better captures the overall band struc- 
ture. It was thus chosen for the further calculations. 
Both parameter sets are reported in Table [H The dif- 
ferences in the parameters is insignificant for the consid- 
erations of the surface superconductivity in the main pa- 
per, and they both give reasonable Fermi velocities (wf,i 
=0.84-10'^ m/s and vp_2 =1.04-10'^ m/s, respectively) . 

DETERMINING THE EFFECTIVE MASS 
DIRECTLY FROM A PARABOLIC FIT TO THE 
DFT BANDS AROUND K 

Additionally, the DFT-calculated dispersion very close 
to the K-point was fitted directly using a parable [Eq. (8) 
in the main text], and a value for the effective mass m* 
and chemical potential were extracted separately for 
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FIG. 4: Fitting of tiie tight-binding parameters for a 20-layer ABC-stacked graphene slab, (a) Comparison of the bands close 
to the Tf-point, as well as the form of the weighing function used in the least squares fit. Inset shows a close-up of the parabolic 
region of the flat band, (b) A comparison of the tight-binding and DFT band structures along the V — K — M lines for the two 
found optima. Black - DFT, blue - tight-binding fit 1, red - tight-binding fit 2. 



both the K — T and K — M directions, in addition to 
a mean value fitted simultaneously to both directions. 
These values were compared to those obtained from the 
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number of layers 



FIG. 5: (a) Effective mass m* obtained from the parabolic 
fit to the DFT bands (black), and calculated from the fit- 
ted tight-binding parameters using [Eq. (8)]. The black open 
triangles refer to m* from DFT correponding to the K — T 
direction, and open squares to K — M direction. Filled black 
circles result from a fit using both directions. The lines are a 
guide to the eye. (b) Tight-binding fits for the parameter a. 
In both plots, blue symbols refer to fit 1 and red symbols to 
fit 2 (Table HJ. 



TABLE I; Tight-binding parameters as a function of the slab 
thickness. Both sets show comparable agreement around K- 
point in the parabolic regime but fit 1 better captures the 
overall shape of the few-layer ABC-stacked graphene slabs. 





fit 1 


fit2 


layers 


70 , 71, 73 , 74 


70 , 71, 73 , 74 


10 


2.62, 0.35, 0.17, 0.06 


3.21, 0.44, 0.16, 0.06 


15 


2.59, 0.35, 0.17, 0.05 


3.21, 0.43, 0.15, 0.05 


20 


2.58, 0.34, 0.17, 0.04 


3.21, 0.43, 0.12, 0.05 



tight-binding fit. Fig.lSJa) shows the m* and /i as a func- 
tion of the slab thickness as well as those calculated from 
Eq. (9) based on the fitted tight-binding parameters 7;. 
The effective mass increases with an increasing number 
of layers, as the breadth of the flat band increases. 

In DFT, there is also a small linear component in the 
parabolic ?A E ~ Ep = l3o{Ak)'^ + /3i(Afc) - fi due to 
the longer-range couplings neglected in the tight-binding 
model. The effect of this component is more important 
further away from K-point, as it decreases when the fit- 
ting region is made narrower. In DFT, the minimum of 
the parabolic bands does not lie exactly at the Fermi en- 
ergy. This finite doping is, however, small, /x/a < 0.1 in 
all DFT calculations. 

The parameter a determines the crossover between the 
BCS- and fiat-band regimes. It is worth noting that even 
though the values for the hopping amplitudes differ in 
the tight-binding parameters sets, both yield comparable 
values for a, as illustrated in Fig. [SJb). Moreover, the 



7 



parameter set given in Ref. |3, 7o = 3.2 eV, 71 = 0.39 eV, 
73 = 0.315 eV and 74 = 0.044 eV, gives a = 0.012 cV, 
showing that this value is quite universal. 
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